LogGabor user guide

Table of content

What is the LogGabor package?

Meanwhile biorthogonal wavelets got a very popular image processing tool, alternative multiresolution transforms have been proposed for solving some of their drawbacks, namely the poor selectivity in orientation and the lack of translation invariance due to the aliasing between subbands. These transforms are generally overcomplete and consequently offer huge degrees of freedom in their design. At the same time their optimization get a challenging task. We proposed here a log-Gabor wavelet transform gathering the excellent mathematical properties of the Gabor functions with a carefully construction to maintain the properties of the filters and to permit exact reconstruction. Two major improvements are proposed: first the highest frequency bands are covered by narrowly localized oriented filters. And second, all the frequency bands including the highest and lowest frequencies are uniformly covered so as exact reconstruction is achieved using the same filters in both the direct and the inverse transforms (which means that the transform is self-invertible). The transform is optimized not only mathematically but it also follows as much as possible the knowledge on the receptive field of the simple cells of the Primary Visual Cortex (V1) of primates and on the statistics of natural images. Compared to the state of the art, the log-Gabor wavelets show excellent behavior in their ability to segregate the image information (e.g. the contrast edges) from incoherent Gaussian noise by hard thresholding and to code the image features through a reduced set of coefficients with large magnitude. Such characteristics make the transform a promising tool for general image processing tasks.

Reference (BibTex Format):

    @article{Fischer07cv,
        Author = {Fischer, Sylvain and Sroubek, Filip and Perrinet, Laurent U. and Redondo, Rafael and Crist{\'o}bal, Gabriel},
        Journal = {Int. Journal of Computional Vision},
        Keywords = {wavelet transforms, log-Gabor filters, oriented high-pass filters, image denoising, visual system},
        Title = {Self-invertible 2{D} log-{G}abor wavelets},
        Url = {http://invibe.net/LaurentPerrinet/Publications/Fischer07cv},
        Year = {2007}}

The log-Gabor transform compared to other multiresolution sche

The log-Gabor transform compared to other multiresolution schemes. a. Schematic contours of the log-Gabor filters implented in Fischer (2007) in the Fourier domain with 5 scales and 8 orientations (only the contours at 78% of the filter maximum are drawn). b. The real part of the corresponding filters is drawn in the spatial domain. The two first scales are drawn at the bottom magnified by a factor of 4 for a better visualization. The different scales are arranged in lines and the orientations in columns. The low-pass filter is drawn in the upper-left part. c. The corresponding imaginary parts of the filters are shown in the same arrangement. Note that the low-pass filter does not have imaginary part. Insets (b) and (c) show the final filters built through all the processes described in Section II. d. In the proposed scheme the elongation of log-Gabor wavelets increases with the number of orientations nt . Here the real parts (left column) and imaginary parts (right column) are drawn for the 3, 4, 6, 8, 10, 12 and 16 orientation schemes. e. As a comparison orthogonal wavelet filters ’Db4’ are shown. Horizontal, vertical and diagonal wavelets are arranged on columns (low-pass on top). f. As a second comparison, steerable pyramid filters [30] are shown. The arrangement over scales and orientations is the same as for the log-Gabor scheme.

Install

Requirements :

  • numpy
  • scipy
  • NeuroTools
  • ipython
  • matplotlib

To install them, use

pip install -r requirements.txt

Install using pip:

pip install LogGabor

Testing filter generation

Tests for LogGabor : filter generation and image filtering

In [1]: 
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
np.set_printoptions(precision=2, suppress=True)
In [2]: 
from NeuroTools.parameters import ParameterSet
pe = ParameterSet('default_param.py')
opts= {'cmap': plt.cm.gray, 'vmin':-1., 'vmax':1., 'interpolation':'nearest', 'origin':'upper'}#, 'figsize':(16,4)}
from LogGabor import LogGabor as LG
In [3]: 
#! this test report is about the LogGabor class:
print LG.__doc__


LogGabor

See http://invibe.net/LaurentPerrinet/Publications/Perrinet11sfn

testing coordinates in Fourier space

In [4]: 
#! let's define the image processing framework
pe.N_X = 32
im = LG.Image(pe)
plt.subplot(111)
#! im.f gives the radial coordinates - note the overlay of coordinates in the matrix
plt.imshow(im.normalize(im.f), **opts)
plt.show()
print im.f[(im.n_x-1)//2 + 1, (im.n_y-1)//2 + 1], im.f[(im.n_x-1)//2 , (im.n_y-1)//2 ]
im = LG.Image(pe)
plt.subplot(111)
plt.imshow(im.normalize(im.f), **opts)
plt.show()
print im.f[(im.n_x-1)//2 + 1, (im.n_y-1)//2 + 1], im.f[(im.n_x-1)//2 , (im.n_y-1)//2 ]

1e-12 1.41421356237


1e-12 1.41421356237
In [5]: 
#! defining input image as Lena
#!-----------------------------
image = plt.imread('database/lena512.png')[:,:,0]
print image.mean(), image.std()

0.499949 0.182783
In [6]: 
pe.N_X = image.shape[0]
im = LG.Image(pe)
image = im.normalize(image, center=True)
print image.mean(), image.std()

-0.000200563 0.403509
In [7]: 
lg = LG.LogGabor(im)
n = image.shape[0]
#! generating some LogGabor filters
#!---------------------------------
#! viewing is done by a dedicated function, show_loggabor :
print LG.LogGabor.show_loggabor.__doc__

None
In [8]: 
#! defining a reference log-gabor 
# (look in the corners!)
sf_0 = image.shape[0]/25. # TODO .1 cycle / pixel (Geisler)
params= {'sf_0':sf_0, 'B_sf': pe.B_sf, 'theta':0., 'B_theta': pe.B_theta}
fig = lg.show_loggabor(0, 0, **params)
In [9]: 
#! translating in the middle
N_X, N_Y = image.shape
fig = lg.show_loggabor(n/2, n/2,  **params)
In [10]: 
#! same params, larger frequency with sf_0
params2 = {'sf_0':sf_0*4, 'B_sf':pe.B_sf, 'theta':0., 'B_theta':pe.B_theta}
fig = lg.show_loggabor(n/4, n/2, **params2)
In [11]: 
#! narrower with B_theta
params3 = {'sf_0':sf_0, 'B_sf':pe.B_sf, 'theta':0., 'B_theta':pe.B_theta/4}
fig = lg.show_loggabor(n/4, 3*n/4, **params3)
In [12]: 
#! broader spectrum with B_sf
params4 = {'sf_0':sf_0, 'B_sf':pe.B_sf*2., 'theta':0., 'B_theta':pe.B_theta}
fig = lg.show_loggabor(3*n/4, n/4, **params4)
In [13]: 
#! a function to explore these parameters:
def lg_explore(param_name, param_range, angle=False):
    for param_ in param_range:
        if angle:
            title = np.str(param_*180/np.pi) + '°'
        else:
            title = np.str(param_)
        print param_name, title
        params_=params.copy()
        params_.update({param_name:param_})
        fig, a1, a2 = lg.show_loggabor(n/2, n/2, **params_)
        #_ = a1.set_xlabel(param_name)
        #_ = a2.set_xlabel(title)
In [14]: 
#! explore these parameters individually:
lg_explore(param_name='phase', param_range=np.linspace(0, np.pi, 6, endpoint=False), angle=True)

phase 0.0°
phase 30.0°
phase 60.0°
phase 90.0°
phase 120.0°
phase 150.0°
In [15]: 
lg_explore(param_name='theta', param_range=np.linspace(0, np.pi, 6, endpoint=False), angle=True)

theta 0.0°
theta 30.0°
theta 60.0°
theta 90.0°
theta 120.0°
theta 150.0°
In [16]: 
lg_explore(param_name='B_theta', param_range=np.pi/np.linspace(4, 14, 6), angle=True)

B_theta 45.0°
B_theta 30.0°
B_theta 22.5°
B_theta 18.0°
B_theta 15.0°
B_theta 12.8571428571°
In [17]: 
#! sf_0  pyramid that we use later in Matching Pursuit
base = 2 #1.618
n_levels = int(np.log(np.max((N_X, N_Y)))/np.log(base))#+1
v_sf_0 = N_X / np.logspace(1, n_levels, n_levels, base=base)
#v_sf_0 =  2. ** np.arange(n_levels)
print v_sf_0, lg.n_x / v_sf_0, len(v_sf_0)==n_levels
lg_explore(param_name='sf_0', param_range=v_sf_0)

[ 256.  128.   64.   32.   16.    8.    4.    2.    1.] [   2.    4.    8.   16.   32.   64.  128.  256.  512.] True
sf_0 256.0
sf_0 128.0
sf_0 64.0
sf_0 32.0
sf_0 16.0
sf_0 8.0
sf_0 4.0
sf_0 2.0
sf_0 1.0

relative distances :

In [18]: 
scaling = np.sqrt(v_sf_0[:, np.newaxis]*v_sf_0[np.newaxis, :])/np.sqrt(N_X)
print v_sf_0
print scaling, scaling.mean()

[ 256.  128.   64.   32.   16.    8.    4.    2.    1.]
[[ 11.31   8.     5.66   4.     2.83   2.     1.41   1.     0.71]
 [  8.     5.66   4.     2.83   2.     1.41   1.     0.71   0.5 ]
 [  5.66   4.     2.83   2.     1.41   1.     0.71   0.5    0.35]
 [  4.     2.83   2.     1.41   1.     0.71   0.5    0.35   0.25]
 [  2.83   2.     1.41   1.     0.71   0.5    0.35   0.25   0.18]
 [  2.     1.41   1.     0.71   0.5    0.35   0.25   0.18   0.12]
 [  1.41   1.     0.71   0.5    0.35   0.25   0.18   0.12   0.09]
 [  1.     0.71   0.5    0.35   0.25   0.18   0.12   0.09   0.06]
 [  0.71   0.5    0.35   0.25   0.18   0.12   0.09   0.06   0.04]] 1.48744418847
In [19]: 
lg_explore(param_name='B_sf', param_range=np.logspace(-1, 0, 5)*pe.B_sf)

B_sf 0.15
B_sf 0.266741911506
B_sf 0.474341649025
B_sf 0.843511987786
B_sf 1.5

Testing on a sample image

In [20]: 
plt.subplot(111)
plt.imshow(im.normalize(image), **opts)
v = plt.axis([0, n, n, 0])
In [21]: 
#! whitening to balance the energy of evey frequency band
white = im.whitening(image)
plt.imshow(im.normalize(white), **opts)
v = plt.axis([0, n, n, 0])
In [22]: 
#! the filtering operation preserves infomation (none is lost...)
plt.imshow(white - im.FTfilter(white, 1.), **opts)
Out[22]:
<matplotlib.image.AxesImage at 0x107632950>
In [23]: 
#! a function to explore these parameters:
def filter_explore(param_name, param_range):
    for param_ in param_range:
        print param_name, param_
        params_=params.copy()
        params_.update({param_name:param_})
        FT_lg = lg.loggabor(0, 0, **params_)
        im_ = im.FTfilter(image, FT_lg, full=True)
        plt.figure(figsize=(18, 6))
        plt.subplot(122)
        im_phi = np.angle(im_)
        plt.imshow(im_phi, cmap=plt.cm.hsv)
        v = plt.axis([0, n, n, 0])
        plt.subplot(121)
        im_a = np.absolute(im_)
        im_a = 2*im.normalize(im_a)-1 # back in the range -1, 1
        plt.imshow(im_a, **opts)
        v = plt.axis([0, n, n, 0])
        plt.show()
#! explore parameters of LogGabors on the filtering:
filter_explore(param_name='theta', param_range=np.linspace(0, np.pi, 6, endpoint=False))

theta 0.0


theta 0.523598775598


theta 1.0471975512


theta 1.57079632679


theta 2.09439510239


theta 2.61799387799
In [24]: 
filter_explore(param_name='B_theta', param_range=np.pi/np.linspace(4, 14, 6))

B_theta 0.785398163397


B_theta 0.523598775598


B_theta 0.392699081699


B_theta 0.314159265359


B_theta 0.261799387799


B_theta 0.224399475256
In [25]: 
print lg.n_x / v_sf_0
filter_explore(param_name='sf_0', param_range=v_sf_0)

[   2.    4.    8.   16.   32.   64.  128.  256.  512.]
sf_0 256.0


sf_0 128.0


sf_0 64.0


sf_0 32.0


sf_0 16.0


sf_0 8.0


sf_0 4.0


sf_0 2.0


sf_0 1.0
In [26]: 
filter_explore(param_name='B_sf', param_range=np.logspace(-1, 1, 5)*pe.B_sf)

B_sf 0.15


B_sf 0.474341649025


B_sf 1.5


B_sf 4.74341649025


B_sf 15.0
In [27]: 
#! TODO: show a tiling of Fourier space with one particular choice of log gabors

#n_levels = 6
##! zooming on some intersing part of the scale space
#for sf_0_ in lg.n_x / 2. / np.logspace(n_levels-3, n_levels-1, n_levels, base=2):
#    print 'sf_0 ' , sf_0_
#    FT_lg = lg.loggabor(0, 0, sf_0=sf_0_, B_sf=pe.B_sf, theta=0., B_theta=B_theta)
#    im_ = im.FTfilter(white, FT_lg)
#    im_ = im.normalize(im_)
#    pylab.imshow(im_, **opts)
#    pylab.show()

Building a pyramid

more book keeping

In [28]: 
%install_ext https://raw.githubusercontent.com/rasbt/python_reference/master/ipython_magic/watermark.py
%load_ext watermark
%watermark

Installed watermark.py. To use it, type:
  %load_ext watermark
30/06/2014 12:39:26

CPython 2.7.7
IPython 2.1.0

compiler   : GCC 4.2.1 Compatible Apple LLVM 5.1 (clang-503.0.40)
system     : Darwin
release    : 13.2.0
machine    : x86_64
processor  : i386
CPU cores  : 8
interpreter: 64bit
In [29]: 
# CSS styling within IPython notebook
from IPython.core.display import HTML
def css_styling():
    styles = open("custom.css", "r").read()
    return HTML(styles)
css_styling()
Out[29]: